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The  work  described  in  this  final  report  has  been  sponsored  by  the  Air  Force  Office  of  Scientific 
Research  through  the  European  Office  of  Aerospace  Research  and  Development  in  London  under 
Grant  AFOSR-91'0118.  The  research  work  period  covered  herein  lasted  from  15.  Nov  1990  to  14. 
Nov  1991  and  presents  the  final  one  after  several  successive  sponsoring  periods.  The  Institut  fur 
Raumfahrtsysteme  (Space  Systems)  of  the  University  Stuttgart  is  therefore  indepted  and  would 
like  to  thank  for  the  longyears  support  on  the  investigation  of  Basic  Processes  of  Plasma  Propul¬ 
sion.  We  also  hope  to  have  contributed  helped  in  understanding  more  profoundly  the  fundamental 
processes  in  and  assessing  more  deeply  the  performance  limits  of  coaxial  magneto  plasmadynamic 
arc  thrusters. 
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Nomenclature 


a 

factor 

A 

cross  section 

S 

magnetic  induction 

e 

electronic  charge 

e. 

specific  energy 

P 

electric  field  strength 

P 

spatial  derivatives 

6 

spatial  derivatives 

h 

Planck’s  constant 

B 

source  vector 

I 

electric  current 

j 

current  density 

J 

Jacobian 

k 

Boltzmann  constant 

m 

mass 

m 

mass  flow  rate 

reduced  mass 

mi^ 

reduced  mass 

n 

particle  density 

n 

normal  vector 

P 

pressure 

P 

pressure  tensor 

9 

flow  variables 

9. 

heat  flux 

Q 

cross  section 

r 

radial  coordinate 

R 

radius 

t 

tangential  vector 

T 

thrust 

t 

time 

T 

temperature 

u 

velocity 

U 

contravariant  velocity 

V 

velocity 

V 

contravariant  velocity 

Vc 

volume 

X 

axial  coordinate 

Xi 

charge  number 

a 

heat  transfer  coefRdent 

plasma  parameter 

^ion 

ionisation  energy 

n 

curvilinear  coordinate 

A 

thermal  conductivity 

M 

viscosity  coefficient 

Mo 

permeability  in  vacuum 

curvilinear  coordinate 

0 

mass  density 

9 

electric  conductivity 

T 

electron  collision  time 

9 

stream  function 

0) 

electron  gyrofrequency 

Subscripts 

€ 

electron 

h 

heavy  particles 

J 

ions 

r 

radial 

ret 

reaction 

s 

surface 

2 

9 

axial 

azimuthal 

V 

paricles 
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1  Introduction 


Since  more  than  one  decade  the  Institute  of  Space  Systems,  University  of  Stuttgart  is  engaged  in 
developing,  testing  and  assessing  the  performance  of  various  types  of  steady  state  coaxial  plasmath- 
rusters.  These  thrusters  [1,  2,  3]  consist  of  a  cathode  rod  in  the  center  and  a  cylindrical  and/or 
nozzle  shaped  discharge  chamber  which  serves  in  the  whole  or  in  part  as  an  anode.  The  gaseous 
propellant  enters  the  arc  discharge  becomes  Joule-  heated  and  ionized  and  assumes  the  lokal  speed 
of  sound  somewhere  in  the  nozzle  throut  or  of  the  constant  cross  sectional  area  of  the  thruster; 
the  plasma  is  than  further  accelerated  eigther  by  termal  expansion  through  a  bell  shaped  nozzle 
and/or  by  axial  electromagnetic  forces.  Depending  on  whether  the  first  or  the  second  acceleration 
mechanisme  is  predominant,  one  speaks  of  an  arc  jet  or  an  magnetoplasmadynamic  (MPD)  thruster 
device.  So  far  only  self  magnetic  thrusters  have  been  investigated  and  developed  at  the  Institut  of 
Space  Systems. 

This  experimental  work  is  accompanied  by  theoretical  investigations[4,  5,  6]  in  order  to  get  a  more 
profound  knowledge  of  the  plasma  behavior  and  the  critical  conditions  which  limit  the  performance 
and  applicability  of  these  thrusters.  Therefore  in  this  report  the  computational  methode  and  the 
more  reasant  code  development  wUl  be  addressed  first.  Afterwards  a  profound  analysis  and  a 
new  explanation  of  the  well  known  performance  limiting  onset  phenomenon  is  presented.  Excellent 
a^eement  is  found  between  calculated  and  measured  onset  conditions  for  the  various  thruster  types 
which  proves  the  validity  of  this  new  explanation. 


2  Code  Development 

Despite  of  the  simplicity  in  design  of  a  coaxial  plasma  thruster  (see  Fig.  1)  the  physics  of  those 
thrusters  is  extremly  complex  and  by  no  means  fully  understood.  For  instance  the  cathode-  and 
anode  attachment  of  the  arc  or  the  arc  boundary  ^ven  by  the  free  stream  transitions  layer  between 
the  electrically  conducting  plasma  and  the  not  conduction  colder  gasstream  should  be  known  better. 
But  despite  of  these  problems  various  codes  have  been  developed  so  far  to  calculate  the  thruster 
performance  based  on  one-dimensional  [1,  7]  ,  quasi  two-dimensional  [4,  8]  and  simplified  two- 
dimensional  codes  assuming  proper  boundary  conditions. 

In  this  report,  a  new  two-dimensional  axisymmetric  code  is  presented,  for  an  MPD  thruster  with 
argon  as  propellant  allowing  ionization  up  to  the  sixth  level.  The  extended  Ohm's  law  is  used  to 
calculate  the  current  contour  lines;  the  electron  temperature  distribution  was  solved  by  including 
the  electron  energy  equation,  whereas  a  two-dimensional  flow  code  is  used  to  obtain  the  velocity, 
pressure  and  heavy  particle  temperature  distributions.  With  this  formulation  the  thermal  nonequi¬ 
librium  in  the  plasma  flow  is  taken  into  account.  However,  the  composition  of  the  plasma  is  cal¬ 
culated  by  means  of  the  electron  temperature  dependent  Saha  equation  for  the  different  ionization 
levels,  the  conservation  of  mass  and  the  electrical  charge  neutrality. 

A  comparison  of  the  numerical  and  experiment^  results  for  different  currents  will  be  presented  and 
discussed  in  this  report. 
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Figure  1:  Principle  Scheme  of  a  Coaxial  Plasmathruster. 


2.1  Outline  of  the  Computational  Method 

The  physical  and  numerical  modeling  of  this  problem  was  developed  at  the  IRS  [4,  7]  For  the 
calculation  of  the  MPD  flow,  a  program  system  is  applied  which  connect  three  different  physical 
fields  to  one  another.  The  three  program  system  parts  are  the  extended  Ohm’s  law  for  plasma 
to  calculate  the  discharge,  the  electron  energy  equation  to  calculate  the  electron  temperature  and 
the  flow  code  to  obtain  the  properties  of  the  flow  field.  The  flow  calculation  was  done  first  by  a 
one-dimensional  flow  code  which  is  still  applied  here  for  the  simulations  of  plasma  thrusters  called 
nozzle  type  MPD  accelerators.  However,  a  fully  two-dimensional  Euler  code  is  used  in  this  report 
for  a  cyUndrical  configuration  of  an  MPD  accelerator,  with  conical  cathode. 

The  three  codes  are  connected  in  the  following  manner:  for  a  given  flow  field  and  an  electron 
temperature  distribution,  the  current,  and  hence  the  magnetic  field  distribution,  was  determined. 
With  the  ohmic  heating,  as  a  result  of  the  current  distribution,  an  electron  energy  calculation  follows 
and  determines  a  new  electron  temperature  distribution.  With  these  results  the  flow  field  equations 
were  integrated.  In  the  next  time  step,  this  new  flow  field  and  the  new  electron  temperature 
distribution  were  taken  to  calculate  the  new  electromagnetic  field  distribution,  and  so  on.  In  the 
following  para^aphes  the  basic  equations  of  the  three  systems  are  presented. 


5 


2.2  Discharge  Code 


The  current  distribution  within  a  selffield  MPD  thruster  is  assumed  to  be  two-dimensional  no 
azimuthal  current  has  to  be  expected.  In  order  to  calculate  the  current  distribution  of  such  an 
arc  discharge,  a  two-dimensional  computer  code  has  been  developed.  The  basic  equation  for  the 
discharge  is  the  extended  Ohm’s  law  for  plasmas: 

j  =  (1) 

Here  is  a  the  electric  conductivity,  We  the  electron  cyclotronfrequency  and  Tc  the  electron  collision 
time.  Rewriting  the  Ohm’s  law  by  means  of  Maxwell’s  equations,  one  obtains  a  vector  equation 
for  the  m^netic  induction  vector  B  in  the  form 

0  =  —(V  X  (- V  X  B))  -  (V  X  («  X  B))  +  —(V  X  i/}{V  x  B)  x  B))  -  V/?Vpe  (2) 

Ho  a  fiQ 

with 

LIT  1 

Dsa  erie 

With  the  Jacobian  J  of  the  transformation  from  cylindrical  to  the  curvilinear  coordinates: 


Bz  dr  ^  dz  dr 

the  matrices  are  formed  in  terms  of  the  derivatives  of  the  cylindrical  coordinates  as  follows: 


c  _ 

’ 


nz  =  -Jz 
rjr  =  J-. 


=  (5) 

where  the  convention  ^^  =  ||,  etc.  is  used. 

The  equation  (2)  yields  with  a  stream  function  ^  =  r  B$  the  elliptical,  partial  differential  equation 
of  2nd  order 

^ 


2al39 

— 3--a/io« 


d9  drj^  ,  _  drj,  ^  ^  drir  ,  _  dtk 
dri  V"  ^  dt} 

r  I  do  do\  /<r»  dpe\  dd  dl3\  2ol39  1 

+  (^' «£  +  ”■  9^)  -  (t  -  j  +  ’^9^ j  +  —  -  H 

r  1  If,  do  do\  (o^  dpe\f,di3  di3\  11 

'■  9« + ’^9,]  *vr--w)  +  ’’•9^)  -  H I 

^  fv  dU  dV\ 
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with  V  and  V  as  the  contravariant  velocities 

U  = 

V  =  +  TjrV  (7) 

The  function  ®(r,  2^)  =  const  now  represents  a  current  contour  line,  since  B  =  3$  is  proportional 
to  where  /(r)  is  the  electric  current  carried  through  the  cross  sectional  area  jrr^.  The  proper 
boundary  conditions  for  follow  from  the  geometry  of  the  thruster  walls  and  electrodes.  At  the 
insulator  inside  the  thruster  and  at  the  inflow  boundary  is  set  to  —  For  the  electrodes  the 
electric  fleld  is  assumed  to  be  normal  to  the  surfaces  •  t  =  0.  In  accordance  with  the  curvilinear 
grid,  where  the  Tj  coordinates  are  perpendicular  to  the  electrodes,  it  is  equal  to  =  0,  and  it 
follows  to: 

^  [^(^r>7r  -6»?2)+  ^  ^  (^r  -  -  Sv)  =  0  (8) 

At  the  other  boundary  sections  and  at  the  symmetry  axis  4”  is  set  to  0.  The  composition  of  the 
Argon  Plasma  is  calculated  as  foUowing. 

Due  to  the  higher  mobility  of  the  electrons,  the  ionization  of  the  argon  plasma  in  MPD  thrusters 
is  dominated  by  the  electron  temperature.  The  Saha  equations  for  the  different  ionization  levels 
[8,  20],  (i=l,2,...,6)  foUow  by 


——Pe  =  -—2 - rr - e  «« 


(9) 


Here  pe  is  <,he  partial  pressure  of  the  electrons,  n,-  the  number  densities,  g,  the  weighting  factor 
taking  into  account  the  degenerated  states  and  €,•  the  ionization  energy  of  the  ionization  level  i. 


The  sum  of  the  partial  pressures  of  all  plasma  components  yields  the  total  pressure 

P  =  (19) 

1/ 

and  the  neutrality  of  charge 

"e  =  (11) 

i 

considering  the  higher  ionization  modes  up  to  the  sixth  ionization  level  the  eqations  9-11  yield  a 
seventh  order  polynom  for  the  electron  density.  In  equation  10  represents  TV  the  electron  temper¬ 
ature  Te  and  the  heavy  particles  temperature,  which  for  aU  heavy  particles  is  equal  T^.  In  eqation 
11  z,  stands  for  the  charge  number. 


The  eqation  for  the  electron  density  is  solved  by  a  Newton  iteration.  With  this  electron  density 
the  other  partial  densities  like  the  argon  atoms  and  the  different  ions  are  determined. 


The  electric  conductivity  of  an  plasma  is  determined  by[9,  21]: 

^  _  3  rW _ e^ne _ 

8  V  2  X<i/(5ie) 


(12) 


Here  nieu  is  the  reduced  mass  and  Te„  the  reduced  temperature.  With  respect  to  the  different 
ionized  levels  the  Gvosdover  cross  sections  follows  by[10,  11,  23]: 


nee^(zf(zi  +  l)) 
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For  Qeu  being  the  cross  section  between  electrons  and  atoms  in  equation  12  the  Ramsauer  cross 
section  is  used  [12,  24].  From  these  eqations  it  is  obviously  that  the  high  ionization  levels  have  a 
strong  effect  on  the  electric  conductivity. 


Te  CkK3 


Figure  2;  Composition  of  an  argon  plasma  as  function  of  electrontemperature  at  1000  Pa. 

For  an  electron  temperature  up  to  80  kK  and  a  constant  pressure  of  1000  Pa  in  Fig.  2,  a  typical 
concentration  distribution  with  respect  to  the  sum  of  heavy  particles  is  presented.  With  the  model 
developed  here  for  higher  ionisation  levels  it  is  obvious  that  calculations  for  typical  MPD  conditions 
are  possible  up  to  an  electron  temperature  of  nearly  80  kK. 


2.3  Electron  Energy  Equation 


The  electron  temperature  has  a  strong  effect  on  the  electrical  and  thermal  conductivity  and  on  the 
electron  density,  which  again  influences  the  discharge  pattern.  Therefore,  a  two-dimensional  code 
for  the  electron  temperature  distribution,  corresponding  to  the  two-dimensional  discharge  code, 
was  written.  The  electron  temperature  distribution  is  determined  by  the  energy  equation  for  the 
electron  component  of  the  plasma. 


V(A,Vre)  -h  "jVTe  +  PeV  .  r  =  -  x;  (T*  -  Th)  -  a 


dt 


(13) 


ret 
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The  subscript  v  represents  the  heavy  particles.  The  first  term  on  the  left  hand  side  represents  the 
heat  flux  in  the  electron  gas,  and  the  second  term  gives  the  convective  change  of  the  electron  gas 
energy  due  to  the  electron  drift.  The  energy  input  due  to  ohmic  heating  is  represented  by  the  first 
term  on  the  right  hand  side.  The  sum  of  losses  due  to  the  energy  transfer  from  the  electron  gas 
to  the  heavy  particles  gas  is  calculated  by  the  second  term  on  the  right  hand  side.  The  reaction 
losses  are  given  by  the  last  term,  where  Cj  is  given  by  Unsold  [8,  20].  The  heat  transfer  coefficient 
Oft,  and  the  thermal  conductivity  Ag  depend  on  the  electron  temperature  and  follow  from  kinetic 
theory[21,  22,  26] 

8\/2  , -  k 

oteu  =  -^Qtuy/i^eukTeu - ; -  (14) 


rrie  +  my 


and  the  thermal  conductivity: 


15  {¥ 

“  8  V  2 


I3i/(=e)  ^vQev  V^v^Ty 


In  the  sum  of  this  conductivity  equation  includes  v  also  the  electron  component.  For  v  =  e  the 
colission  cross  section  Qee  can  be  set  equal  to  the  Gvosdover  cross  sectionflO,  23].  If  Qey  being  the 
cross  section  between  the  electrons  and  ions,  the  Gvosdover  cross  section  is  used  with  respect  to 
the  different  ionization  levels,  and  for  Q^y  being  the  cross  section  between  electrons  and  atoms  the 
Ramsauer  cross  section  is  used  [12,  24]. 

With  respect  to  the  rotational  symmetry,  equation  13  results  in  the  following  elliptical,  partial 
differential  equation  of  2nd  order: 

._2  ,  .  /,2  .  ,2x^'rg  ,  _  .  ,  _  ,  . 

(^r  "f"  ‘Hz)  Qfj2  ^ 


P 

4AggTgV2 


E/rr  'T  \  I  J 

^  n,nya,y  {Th  -  Tg)  +  -  -  e.-  j  - 

3ngfc  (  dTe  aTg  ^  aTg  dT,\ 

—  (^ur,—  -  ur,—  +  vz,—  -  vz,—)  - 

n,kT,  .  ^  ,]  aTg  aTg  ^ 

—J—  ® 


Here  ^  and  rj  indicate  the  partial  differentiation  to  the  curvilinear  and  rj—  directions.  J  is  the 
Jacobian  of  the  transformation  and  Oi  ,  02  are  metric  terms. 

At  the  outflow  boundary  Tg  is  set  to  a  value  in  accordance  with  measurements[13,  25].  The  solid 
bodies  of  the  thruster  are  treated  as  thermal  insulators.  Therefore  the  boudary  condition  inside 
the  thruster  is  given  by  VTg  •  n  =  0,  where  n  is  the  normal  vector  of  the  surfaces.  Due  to  the  axial 
symmetry,  ^  =  0  on  the  axis.  At  the  inflow  boundary  Tg  is  set  to  a  constant  value  of  7000  K  in 
accordance  with  the  measurements. 


2.4  Flow  Field  Code 


Since  the  flow  field  computation  requires  the  most  calculation  time,  a  one-dimensional  flow  code  is 
developed[14]  for  numerous  geometrical  and  paramater  variations.  In  case  of  the  one-dimensional 
flow  field  modeling,  the  equation  of  state  together  with  the  continuity  and  Bernoulli  equations  yield 
the  following  expansion  equations: 
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-  inside  the  thruster  (isothermal  expansion); 

-  outside  the  thruster  (adiabatic  expansion): 


(17) 


(18) 


where  v^o,  Peo  a^nd  Peo  ‘•■re  the  velocity,  pressure  and  density  of  the  plasma  at  the  thruster  nozzle  exit. 


For  the  description  of  the  two-dimensional,  axisymmetric  flow,  the  following  nonlinear  hyperbolic 
system  of  differential  equations  with  curvilinear,  cylindrical  coordinates  is  used; 


(19) 


The  indices  t,  f/  indicate  the  partial  differentiation  with  respect  to  time  and  to  the  and  p- 
directions.  The  first  three  terms  are  used  in  the  usual  fluid  dynamic  manner,  where  q  is  the  flow 
variables  vector  and  F  and  G  are  the  spatial  derivatives  vectors  [15,  27].  The  source  vector  H  is 
developed  to: 

pvfr 


^(91  = 


puv/r  + 
pv^lr  +  ^  {m 


2W  +  (f  +  f) 

du  ,dV\ 

dr,  +'5?’) 


2f+( 


}  +  jriBe 

]+kBe 


L  (p  -I-  e)v/r  +  Vq,-Ylu  n^neOe^iTh  -  T^) 

The  first  terms  in  the  source  vector  transform  the  plane  two-dimensional-  into  a  cylindrical  calcu¬ 
lation.  The  first  additional  terms  in  the  impulse  equations  represent  the  stress  tensor  of  the  plasma 
flow,  where  the  viscosity  coefficient  p  is  given  by  [9]: 

minikTi 


(20) 


Since  the  temperature  of  the  heavy  particles  are  equal,  it  is  T,  =  Ti„  =  Tfc  and  m/„  is  the  reduced 
mass  lor  /,  i/  as  heavy  particles  subscripts.  The  last  terms  in  the  impulse  equations  represent  the 
j  X  B  forces.  In  the  energy  equation  the  second  source  vector  term  q^  represents  the  heat  flux 
vector  due  to  the  heavy  particles  of  the  MPD  flow,  which  is  given  by: 


9.  =  E-M“Vn,-hVT,) 


(21) 


The  heat  conductivity  coefficients  are  given  by: 

A  -  —  (22) 

'  8  V  2  !:.(=/)  ^  ’ 

The  transfer  of  Joule’s  heat  is  contmned  in  the  temperature  compensation  between  the  electron 
and  the  heavy  particles  temperature,  which  is  represented  by  X;n„neOe^(Te  -T^)  in  the  equation  of 
energy  conservation.  The  transport  coefficients  are  derived  in  a  similar  manner  as  in  [16,  28]  and 
more  detailed  derivations  for  the  coefficients  are  described  in  [9,  21].  The  cross  sections  Qi^  for  all 
coefficients  are  taken  from  the  references  [10,  12,  23,  24]. 
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2.5  Discussion  of  Results 
2.5.1  Nozzle-type  MPD  Thrusters 

A  main  goal  for  developing  tne  code  with  a  one-dimensional  flow  calculation  was  to  obtain  a  design 
instrument  which  provides  a  fast  insight  into  the  current  distribution  within  a  continuously  running 
MPD  thruster.  The  next  figures  show  the  distributions  of  the  main  parameters  calculated  for  three 
nozzle-type  throat  geometries  (see  Fig.  3).  The  DT2  geometry  has  the  narrowest  throat  diameter 
of  24  mm,  whereas  the  DT5  has  a  bigger  throat  diameter  of  30  mm,  and  the  DT6  thruster  doesn’t 
have  a  throat  at  all,  but  has  a  cylindrical  part  with  a  diameter  of  36  mm  from  the  cathode  up 
to  the  beginning  of  the  nozzle.  With  these  three  geometries  it  can  be  expected  that  at  the  DT2 
thruster  the  highest  values  of  current  density,  and  thus  the  highest  temperatures,  occur  in  the 
nozzle  throat.  According  to  the  experiences  in  [14,  17,  18,  26]  in  the  DT2  thruster  also  the  highest 
voltage  is  predicted. 


With  quasi-steady  pulsed  thrusters  it  is  possible  to  measure  the  current  contour  lines,  but  because 
of  the  high  heat  load  of  the  probes  it  is  not  yet  possible  to  compare  these  with  experimental  data 
on  the  continuously  running  thruster.  This  means  there  is  a  need  for  a  calculation  code  in  order 
to  optimize  steady  state  MPD  thrusters.  In  Fig.  4  the  calculated  current  contour  lines  are  shown 
for  the  different  throat  geometries  at  3000  A  with  a  mass  flow  of  0.8  g/s  argon. 
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Figure  4:  Current  contour  lines  of  the  DT-geometries. 

The  electron  temperature  distribution  within  and  outside  the  thruster  at  a  current  level  of  2  kA 
is  shown  in  Fig.  5.  Against  any  expectation,  the  highest  electron  temperature  is  calculating  for 
the  DT6  thruster.  This  is  due  to  the  fact  that  the  calculation  which  is  necessary  for  the  coupling 
of  the  electron  and  heavy  particles,  as  described  in  equation  13,  has  not  been  considered  due  to 
numerical  difficulties. 
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Figure  5:  Electron  temperature  of  the  DT-geometries. 

The  potential  lines  within  the  discharge  at  3  kA  are  plotted  in  Fig.  6.  An  integration  across  these 
lines  yields  the  discharge  voltage.  As  expected,  the  DT2  thruster  has  the  highest  and  the  DT6 
thruster  the  lowest  voltage  value.  This  calculated  discharge  voltage  is  compared  with  the  measured 
voltage  in  the  current  range  between  1  and  5  kA  in  Fig.  8. 
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Figure  6:  Potential  lines  of  the  DT-geometries. 

Fig.  7  shows  the  ionization  distribution  of  the  argon  plasma  up  to  to  the  5th  ionization  level.  This 
ionization  number  shown  here  is  calculated  like  the  electron  concentration  in  Fig.  2,  with  respect 
to  the  sum  of  the  heavy  particles  to  a  =  nd  n„.  In  comparison  to  the  electron  temperature  one 
receives  too  high  ionization  values  in  the  DT6  thruster  according  to  the  lacking  coupling  between 
the  energy  of  the  electrons  and  heavy  particles. 
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Figare  7:  Ionization  distribution  of  the  DT-geometries. 

The  measured  discharge  voltages  include,  in  contrast  to  the  calculation,  the  electrode  faU  voltages. 
In  the  case  of  the  cathode  this  voltage  drop  equals  about  the  work  function  of  the  cathode  material, 
which  means  with  of  thoriated  tungsten  between  1  and  3  volts,  and  is  nearly  independent  of  the 
discharge  current.  The  anode  fall  voltage  depends  strongly  on  the  current;  starting  with  several 
negative  volts  at  currents  up  to  about  3  kA,  the  anode  voltage  rises  steeply  in  the  onset  (unstable 
mod)  reffon  [1,  7]. 
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Figure  8:  Calculated  and  measured  discharge  voltage. 


According  to  the  flow  conditions  and  the  known  electromagnetic  force  configuration,  the  thrust  can 
be  calculated.  The  thrust  of  an  MPD  thruster  is  the  sum  of  all  gas  dynamic  surface  forces  and  the 
electrom^netic  volume  forces.  Hence  it  is 

f  =  +  Pj -dA- j^Jx  B  dV  (23) 

where  A,  represents  the  surface  of  all  internal  walls  and  V  is  the  current  carrying  volume. 


Figure  9:  Calculated  and  measured  thrust. 

The  calculated  thrust  is  compared  with  experimental  data  in  Fig.  9.  For  all  thrusters  with  the 
three  different  geometries  the  calculated  thermal  thrust  was  verified  with  the  thermal  efficiency 
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from  the  experiment.  The  data  shows  a  very  good  agreement  spedaly  for  the  DT2  thruster  for  the 
complete  presented  current  range. 

The  results  for  the  DT2  and  DT5  thrusters  show  a  good  correspondence  with  experimental  values. 
The  calculated  results  for  the  DT6  thruster  show  a  bigger  discrepancy  to  measured  values  because 
of  the  lacking  coupeling  between  the  electron  and  heavy  particles  energy.  For  this  geometry  a 
better  correspondence  to  experimental  values  can  be  achieved  by  improving  the  one-dimensional 
flow  model.  Another  possibility  would  be  a  two-dimensional  flow  modelling  containing  a  comlete 
coupeling,  like  used  for  the  ZT3  thruster,  as  presented  in  the  next  chapter,  but  then  with  a  major 
increase  of  computation  time. 


2.5.2  Cylindricsd  MPD  Thruster 

The  next  figures  show  results  for  a  cylindrical  MPD  thruster  geometry,  the  ZT3-IRS  (see  Fig.  10). 
Here  the  flow  field  was  calculated  by  the  Euler  equations. 

The  flow  inlet  boundary  conditions  were  iterated  to  coincide  with  an  experimentally  obtained  cold 
gas  thrust  of  1  iV  at  a  mass  flow  rate  of  2gls.  For  a  ^ven  current  of  6  kA  the  computation  yields  the 
current  denuty  distribution  as  illustrated  in  Fig.  11.  The  calculated  current  contour  distribution 
corresponds  with  the  current  contour  distribution  at  a  continuing  mode  if  the  cathode  is  hot  glowing 
and  is  emitting  electrons  thermionically  along  its  complete  length  (see  Fig.ll). 


Figure  10:  Cylindrical  MPD  thruster  ZT3-IRS. 
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Figure  13:  Heavy  particles  (top)  and  electron  (bottom)  temperature  contours  for  10  kA. 


The  heavy  particle  (top)  and  the  electron  temperature  distribution  (bottom)  within  and  outside  the 
ZT3  thruster  is  shown  in  Fig.  12  for  a  pven  current  of  6  kA.  The  maximum  temperature  value  of 
the  heavy  particles  occurs  at  the  tip  of  the  cathode  on  the  symmetry  axis.  The  majdmum  electron 
temperature  value  occurs  also  at  the  tip  of  the  cathode,  but  there  is  also  a  week  increase  of  the 
electron  temperature  at  the  beginning  of  the  anode. 

In  Fig.  13  the  heavy  particle  (top)  and  the  electron  temperature  distribution  (bottom)  is  shown  for 
a  ^ven  current  of  10  kA.  The  maximum  temperature  value  of  the  heavy  particles  and  the  electrons 
is  also  at  the  tip  of  the  cathode,  but  there  is  now  a  remarkable  increase  of  the  electron  temperature 
at  the  beginning  of  the  anode. 

The  pressure  map  (Fig.  14)  demonstrates  the  expansion  flow  with  the  flow  described  in  section 
2.3,  where  the  relatively  low  pinch  effect  was  taken  into  account. 
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Figure  14:  Pressure  distribution. 


The  corresponding  velocity  distribution  to  the  temperatures  and  density  distribution  is  shown  in 
Fig.  15  as  a  vector  graph.  It  shows  a  high  increase  from  the  inflow  boundary  downstream.  The 
radial  velocity  components  inside  the  channel  are  quite  small,  which  is  also  an  indication  for  a 
relatively  low  pinch  effect  there. 
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Figure  15:  Velocity  vector  distribution. 

The  calculated  thrust  is  5.2  N  for  4  kA,  8.1  N  for  6  kA,  11.7  N  for  8  kA,  16.4  N  for  10  kA  and 
22.1  N  for  12  kA.  The  measurements  for  the  ZT3  thruster  had  not  been  finished  by  the  time  the 
report  was  completed. 
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2.6  Conclusions  on  Code  Development 


The  codes  described  in  this  report  were  developed  and  tested  in  order  to  gain  a  more  profound 
understanding  of  the  fundamental  processes  occuring  in  MPD  thrusters  to  predict  thruster  perfor¬ 
mances  under  various  operating  conditions  and  to  end  up  with  reliable  design  criteria. 

The  presented  numerical  analysis  of  an  nozzle  type  MPD  thruster  was  applied  to  rather  com¬ 
plex  geometries.  Despite  the  shortcomingsof  this  code,  namely  the  decoupled  computation  of  the 
electromagnetic  and  flow  equations,  the  one-dimensional  flow  model  and  the  numerical  difficulties 
concerning  the  energy-coupling  between  electrons  and  ions  especially  for  thrusters  with  a  large 
nozzle  throat,  it  yields  valuable  results: 

•  Calculated  integral  values,  such  as  thrust  and  discharge  voltage  fit  the  measured  ones  very 
well. 

•  The  current  patterns  of  the  discharge  coincide  with  the  measured  data  fairly  well. 

Considering  the  saving  in  computation  time  this  code  is  a  useful  design  instrument. 

The  two-dimensional  flow  field  calculation  code  with  the  curvilinear  koordinates  permits  an  im¬ 
proved  simulation  of  geometry  dependencies.  This  code  was  validated  for  the  cylindrical  thruster 
ZT3-IRS  and  showed  good  agreement  with  the  experiments  (Tab.  1)  for  the  voltage  and  the  cur¬ 
rent  distribution  on  the  three  anode  segments.  Therefore  in  a  later  research  period  also  the  nozzle 
type  thruster  designs  will  be  investigated  with  this  two-dimensional  curvilinear  code. 


S  Thruster  Stability  Investigation 


One  of  the  most  critical  limiting  effects  of  the  self-magnetic  MPD-thrusters  turns  out  to  be  the  so 
called  ‘‘Onset  Phenomenon”  [17, 18, 19, 20, 21]  a  plasma  instability  which  causes  severe  degradation 
of  the  thruster  performance  at  high  power  levels.  According  to  experimentals  done  with  the  steady 
state  nozzle  type  thrusters  (see  Fig.  3)  for  a  given  flow  rate  the  dc  voltage  rises  steeply  if  one 
exceeds  a  critical  current  (see  Fig.  16);  simultanously  voltage  fluctuations  can  be  observed  the 
intensity  of  which  grow  rapidly  by  further  increasing  the  current.  This  effect  finally  leads  to  a 
strong  degradation  of  the  thruster  by  locally  overheading  the  anode  and  irreversibly  damaging  the 
device.  Similar  observations  could  be  observed  on  pulsed  thruster  devices  [18, 21,  24,  25].  However, 
for  a  pulsed  thruster  experiment  a  correct  measurement  of  the  effective  mass  flow  is  extremly 
difficult  since  in  fact  a  quasi  steady  condition  near  the  cathode  attachment  and  especially  of  the 
plasma  cold  gas  interface  within  the  discharge  chamber  remains  highly  questionable. 


Current  [  kA  ] 


Figure  16:  Voltage  vs.  current  curve  for  a  mass  flow  rate  of  0.8  g/s  of  argon 

Several  different  explanations  [7,  8,  9,  10,  12,  13,  22,  23,  26,  27,  28,  29]  have  been  put  foreward  so 
far  and  for  most  of  them  their  analytical  results  agree  fairly  well  with  experimental  observations. 
In  this  report  a  new  explanation  of  the  ‘Onset- Phenomenon’  is  presented  which  is  based  on  a  ‘run 
away  Joule  heating’  mechanism  of  the  current  carrying  plasma  channel.  This  explanation  is  in 
excellent  agreement  with  the  experimentally  found  onset  conditions  measured  on  our  continously 
running  nozzle  type  thruster  devices  and  therefore  worth  to  be  mentirned. 

It  is  based  on  the  ‘WgTe-dependency’  of  the  electron  heatflux  (w*  =  cyclotron  frequency  of  electrons, 
Tg  =  average  collision  time)  and  the  fact  that  with  increasing  a  contraction  and  overheating  of 
the  discharge  channel  may  occur.  Such  a  process  consequently  would  lead  to  a  ‘run  away  heating’ 
effect  and  as  shown  elsewhere  [14,  30]  to  an  instability  of  the  rotationally  symmetrical  channel 
conjuration.  Based  on  this  stability  theory  the  new  quasi  stable  configuration  of  the  constricted 
discharge  channel  assumes  a  helical  shape. 
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3.1  Run  Away  Heating  Process 


Within  the  constrictor  or  nozzle  throat  of  a  self  magnetic  MPD-thruster  the  plasma  conditions  are 
characterized  by  an  axial  flow  at  Mach  1  and  the  fact  that  the  Joule  heat  input,  due  to  the  axial 
electric  current  which  primarily  heats  up  the  electron  gas  is  balanced  by  the  radial  conduction  heat 
flux  of  the  electrons  and  by  the  energy  transfer  to  the  heavy  particle  components.  Since  the  pressure 
in  those  devices  is  well  below  1  bar,  one  can  neglect  radiation  losses  and  moreover  may  assume  that 
the  rotationally  symmetric  discharge  column  fills  up  the  entire  constrictor.  If  the  pressure  or  the 
mass  flow  would  become  too  high  (arcjet  conditions)  the  cold  gas  region  between  discharge  column 
and  constrictor  wall,  and  hence  the  axial  convection  effects  in  the  transition  regime  between  outer 
arc  column  and  cold  gas  should  be  accounted  for.  The  governing  equation  of  the  constrictor  follows 
now  by  the  simpUfied  energie  equation  of  the  electrons  (see  eq.(13))  as 


Now,  within  the  discharge  the  conduction  heat  flux  of  the  heavy  particle  components  is  approxi¬ 
mately  a  factor  =  3, 7  •  10“^  smaller  than  that  of  the  electrons.  One  can  therefore  propose 
that  within  the  downstream  area  of  the  constrictor  the  heavy  particles  assume  the  temperature  of 
the  directly  Joule  heated  electrons.  Consequently  in  the  downstream  area  of  the  constrictor  the 
second  term  on  the  right  site  of  equation  (24)  may  be  neglected.  The  solution  of  the  remaining 
equation  yields  the  following  electron  temperature  profile  [5]. 


I  «(« 
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where 


Te  =  1.513  •  10Vo7  j  (26) 

is  the  maximum  electron  temperature  at  the  center  line  degree  Kelvin  if  is  given  in  .  A 
stands  for  the  electric  current  which  is  transported  through  the  plasma  column  of  radius  Tc-  This 
column  radius  is  taken  to  be  identical  to  the  constrictor  radius.  The  temperature  profile  has  been 
derived  as  a  solution  of  the  electron-energy  equation,  (eq.  (24)  with  Te  =  T  ),  by  modeling  the 
current  density  j  through  a  paraboloid  of  grade  n. 


where  the  maximum  current  density  follows  by  the  relation 


n  -f-  2  / 

j  _  _  __ 


23 


The  product  of  the  two  factors  Jo  and  7  in  eq.(26)  is  a  dimensionless  number  between  about  one  and 
three.  Herein  the  factor  Jq  turns  out  to  be  a  function  of  grade  n  of  the  current  density  paraboloid. 
In  Table  1  Jq  is  given  for  various  values  of  n. 


Table  1:  Jq  (n) 


n= 

1.5 

2 

3 

4 

10 

100 

1.297 

1.246 

1.185 

1.149 

1.070 

1.008 

The  dimensionless  factor  Jq 

The  factor  7  depends  on  the  type  of  gas  and  is  dimensionless,  slowly  varying  function  of  fg  and 
the  pressure  p,  or  of  ^  and  p  within  the  constrictor  for  an  argon  plasma.  This  factor,  7,  can  be 
presented  by  [5] 


lA  ^ 


0.853 

~  — 
1  -  0.00189  [1  -  0.22  •  In  (j^)]  Jq  " 


(29) 


Here  ^  must  be  taken  in  and  p  in  [Pa].  Equation  (29)  is  valid  within  a  pressure  range  of 

10^ Pa  <  p  <  10® Pa  and  an  -range  of  5  •  10“*^  <  ;^  <  5  •  10®^,  provided  the  maximum 

temperature,  fg  stays  within  the  lower  and  upper  limit  of  15000  K  and  40000  K,  respectively.  The 
tolerated  error  within  these  ranges  is  then  less  than  ±  5%. 


Figure  17:  Temperature  Profile  for  Different  Grades  n  of  the  Current  Density  Paraboloides  within 
the  Arc  Core 
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The  dimensionless  temperature  profile  according  to  eq.(25)  is  now  plotted  for  various  current  den¬ 
sities  or  ‘grades  n’  in  Fig.  17.  As  can  be  seen  the  temperature  profile  varies  only  slightly  with  n, 
a  fact  that  justifies  the  method  of  solving  the  constrictor  condition,  i.e.  by  modeling  the  current 
density  across  the  constrictor  area. 

So  far  the  derived  temperature  profile  was  calculated  on  the  assumption  that  no  magnetic  field 
effects  are  present  or  of  importance.  As  it  is  known,  however,  within  and  around  the  current 
carrying  arc  column  a  self  magnetic  field  or  induction  field  B  is  inducted  pointing  in  the  azimuthal 
ot  6  -  direction  (see  Fig.  18).  Based  on  the  current  density  profile  of  eq.(27)  with  eq.(28)  the  radial 
depensity  of  this  azimuthal  magnetic  induction  field  is  given  by 


Figure  18;  Illustration  of  the  Self  Magnetic  Field  within  a  Plasma  Thruster 


As  a  consequence  of  this  azimuthal  magnetic  field  which  is  normal  with  respect  to  the  radial 
temperature  and  density  gradients,  the  radial  heat  flux  especially  that  due  to  electrons  will  be 
affected.  Based  on  a  rigorous  gas  kinetic  perturbation  approach  which  takes  into  account  such  a 
megnetic  field  [5]  ,  it  can  be  shown  that  the  heat  conduction  coefficient  as  given  by  eq.(15)  must 
be  corrected  by  a  tensorial  factor  which  depends  on  the  product  of  the  cyclotron  frequency  Wg  and 
the  average  collision  time  of  an  electron  r^.  Now  it  is 


(31) 


directly  proportional  to  the  magnetic  induction  field  B,  and  the  average  collision  time  as  gas  kinetic 
expression  given  by  [5] 
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(32) 


T.  =  - 
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Ti^ 
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Within  a  magnetic  field  the  heat  conduction  coefficient  therefore  becomes  a  tensor  which  follows 

by 
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(33) 


where  now  ^2  =  0,  By  =  —B,  B^  =  0  and  u)i  =  =  —  ^  for  i  =  e.  Since  within  the  constrictor 

only  a  radial  temperature  gradient  -  i.e.  parallel  to  the  x-direction  (see  Fig.  18)  -  is  present  the 
conduction  heat  flux  vector,  of  the  electrons  now  becomes 


9e  — 


{ 


StnUeTe  ^ 
- «r 


1  -  COSUeTe 
WeTe 


(34) 


where  u,  and  Ux  are  the  unit  vectors  in  the  r-  and  z-  direction,  respectively.  If  the  self  magnetic 
field  {B  (JeTe  ^  0)  is  accounted  for,  the  heat  flux  vector  not  only  has  a  radial  component  but  also 
an  axial  one.  However,  this  axial  heat  flux  does  not  lead  to  a  net  heat  loss  within  the  constrictor 
since,  ‘^what  flows  axially  out  flows  axially  in”.  Yet,  the  radial  heat  flux  would  lead  to  a  net  heat 
loss  if  it  were  not  steadily  compensated  by  Joule  heating.  This  radial  head  flux  now  depends  on  the 
magnetic  field  and/or  on  the  product  of  cyclotron  frequency  and  the  average  collision  time,  UeTe, 
according  to  the  factor 


sinUeTe  ^  _  (UeTef  {(^eTe)* 

UxTe  6  120  '  ■  ' 


(35) 


For  very  small  WeTe  -values  (  We^'e  <  0.1  )  this  factor  differs  only  by  a  negligible  amount  (-0.002) 
from  one.  The  magnetic  field  eflfect  can  therefore  be  neglected  as  long  as  WeTe  <  0.1.  If  however  WeTg 
exceeds  a  critical  value,  (weTe)cr,(  0-5  •  •  -l-O),  this  factor  eq.(35)  becomes  noticeably  smaller  than 
one  and  the  magnetic  field  effect  on  the  radial  heat  conduction  of  the  electrons  must  be  accounted 
for.  As  one  may  realize,  a  noticeable  reduction  of  the  radial  heat  conduction  consequently  leads  to 
a  temperature  increase  within  the  core  of  the  discharge  channel. 

Since  now  the  self  magnetic  field  B(r)  and  the  electron  cyclotron  frequency  Ug  have  a  maximum  at 
a  given  radius  tm,  according  to  (see  eq.(30))  the  optimization  rule 


the  electron  experimence  the  strongest  megnetic  field  and  have  the  highest  cyclotron  frequency,  Wg 
at  the  radius  of 


rM  = 


(36) 
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Therefore  inside  a  radius  of  about  the  heat  conduction  losses  will  decrease  at  first  and  simul¬ 
taneously  increase  the  temperature.  However,  with  increasing  temperature  the  collision  frequency 
of  the  electrons  (see  eq.(14))  rises  to  such  an  extent  that  UeTe  increases  even  further  and  the  con¬ 
duction  heat  loss  of  the  electrons  are  almost  turned  off.  The  consequence  is  a  run-away  heating 
within  the  channel  core  (r  <  rjvf )  until  a  new  balanced  condition  is  achieved.  This  new  condition 
may  be  characterized  by  the  fact  that  the  radial  conduction  loss  due  to  electrons  ceases  (  o^e^e  =  ^ 
)  or  is  negligible  and  that  the  Joule  heat  -put  into  the  electrons-  is  balanced  by  the  energy  transfer 
to  the  heavy  particles.  This  energy  is  then  compensated  by  the  heat  conduction  loss  of  the  heavy 
particles.  The  new  temperature  profile  should  therefore  be  determinated  by  the  much  smaller  heat 
conduction  losses  of  the  heavy  particles  than  by  the  classical  electron  conduction  losses.  Neverthe¬ 
less,  based  on  this  discussion  it  may  be  concluded  that  if  UeTe  taken  at  a  radius  of  exceeds  a 
critical  value  of  (i^eTe)crit  ^  0.5  •••1.0,  the  temperature  profile  should  increase  in  the  core  region 
(r  <  rjvf)  and  decrease  outside  (r  >  r^),  causing  a  contraction  of  the  discharge  cross  section  (see 
Fig.  19). 


Figure  19:  Illustration  of  the  (jj^Te  -effect  on  the  temperature  profile 


3.2  Calculation  of  (wgTe),^ 

One  now  obtains  according  to  eq.(31)  and  eq.(30)  the  cyclotron  frequency  &i  r  =  tm  (see  eq.(36)) 
by 


where 


(  ^  ^  ^  T 

—  ^eAfax  —  o  ^ 
TTlg  ZK  Tq 


(37) 


Jz  = 


(38) 


is  a  function  of  the  grade  n  of  the  current-density  paraboloid  (eq.(27)).  J3  varies  between  1,1  and 
1,0  if  n  changes  from  1,5  to  100,  respectively. 

The  collision  time  according  to  eq.(32)  can  now  be  rewriten  in  the  form 
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where  Qee  follows  by  the  Gvosdover  cross  section  with  =  —1.  In  the  denominator 

of  eq.(39)  the  sum  with  the  terms  proportional  to  next  to  one  can  now  be  neglected.  ( 

rrieT^  <  m^Te ;  <  1  as  long  as  one  does  not  consider  plasmas  of  a  low  degree  of  ionization, 

i.e.  low  temperatures  ;  lOOOOK).  Together  with  eq.()  one  therefore  obtains 


Te 


3 

^r2 


By  means  of  the  definition  of  the  degree  of  ionization 


and  the  pressure  relation 


a  = 


We 


(40) 


(41) 


P  =  53  +  ‘f^’ekTe 

*'(/«) 


(42) 


the  electron  density  can  be  expressed  by 


a  p 
“  1  +  0^ 


(43) 


if  constrictor  conditions  with  Tg  ^  T  are  assumed.  On  replacing  the  electron  density  rig  in  eq.(40) 
by  eq.(43),  the  average  collision  time  follows  as  a  function  of  Tg  and  p. 


_ 24\/2xt _  1  +  Q  epy/m^  {kTg)i 

3.80  •  10-**  1  +  a  Te^ 

/n|5.705  10-3(l^)*r,|  "  ^ 


(44) 

(45) 


In  order  to  obtain  the  collision  time  at  the  effective  channel  radius  tm,  one  has  to  take  for  the 
electron  temperature  in  eq.(45)  Tg  =  Tg  (tm).  This  temperatrue  follows  according  to  eq.(25)  and 
(36)  by 


Tgiru) 


=  n|i- 


(n  4-  2)(n  -f  3) 
(n  +  l)(n  +  4) 


=  fgJ2 


(46) 
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where  J3  is  a  dimensionless  factor  of  the  order  0,6  which  depends  on  the  grade  n.  Inserted  into 
eq.(45)  one  obtains  the  effective  collision  time  by 


(a;ereU  =  3.76.1O-«/^(0'  (47) 

Herein  the  p  must  be  taken  in  Pascal,  I  in  Ampere  and  Tc  in  meter  in  order  to  obtain  (u;eT’e),jj/ 
dimensionless.  The  factor  /  is  of  the  order  one,  depending  on  the  type  of  gas  and  is  a  slowly  varying 
function  of  n,p  and  Tg  or  of  n,  p  and  It  is 


_ J3  _ 1  +  0 

in  (5.705  “ 


(48) 


At  lower  temperatures,  Tg  {tm)  =  J^Te,  and  higher  pressures  the  degree  of  ionization  may  become 
fairly  small  (a  <  1)  such  that  /  would  become  quite  large.  This  range,  however,  is  not  of  interest 
here  and  would  result  in  wrong  /  -values  since  the  electron-atom  collision  were  neglected  in  Tg  and 
the  transport  coefficients.  The  formula  of  (48)  should  therefore  only  be  used  for  a  >  0.6. 

It  should  be  remarked  here  that  eq.(47)  has  been  derived  under  the  assumption  that  the  radial  heat 
conduction  is  still  not  noticeably  affected  by  the  UgTg  -effects,  i.e.  the  factor  with  which  the 

radial  conduction  heat  flux  should  be  multiplied  was  simply  assumed  to  be  approximately  equal  to 
one  (see  previous  section).  Now  for  (we^e)rM  ~  the  relative  error  in  the  heat  flux  calculation 
amounts  to  only  1%  and  even  for  {ugTg)^f^  =  0.75  it  is  less  than  10%.  It  may  therefore  be  concluded 
that  if  {ugTg)^f^  stays  somewhat  below  about  0.75,  eq.(47)  can  be  applied  to  determine  {ugTg)^!^. 
On  the  other  hand  the  {ugTg\j^  -value  calculated  according  to  eq.  (47)  should  never  exceed  a 
critical  value,  (wgr*)^,-.  (<  0.75)  in  order  to  avoid  run  away  heating  and  sudden  temperature  and 
current  density  concentration  towards  the  center  fine  of  the  constrictor.  Instead  of  relating  (weTe)^^ 
to  the  pressure  p  as  given  in  eq.(47),  p  may  be  replaced  by  the  mass  flow  rate,  rh  of  the  thruster. 
Provided  that  the  discharge  fills  the  entire  cross  section  of  the  constrictor  (xr^)  the  mass  flow  rate 
is  obtained  by 


rhc  =  2x  /  pvrdr  (49) 

Jo 

where  p  is  the  local  density  given  by  the  relation  (see  also  eq.(42)) 

where  Mo  is  the  atomic  weight  of  the  propellant,  K  is  the  universal  gas  constant  and  T  (=  Tg)  is  the 
temperature  of  the  heavy  particles  and  within  the  constrictor  set  equal  to  the  electron  temperature. 
The  velocity  v  can  be  replaced  by  the  local  speed  of  sound  under  isothermal  conditions  and  follows 
as 


(50) 


V 


(51) 
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Hence,  the  axial  mass  flow  through  the  constrictor  can  be  calciilated  by 


=  2Tr-=J _ rdr 


Neglecting  the  magnetic  pinch  effect  -  which  is  indeed  allowed  here  -  one  obtains  with  n=4  for  the 
mass  flow  rate 


The  =  TTfe 


I  Mo  1.15 

nfVTT& 


where  d  is  an  average  degree  of  ionization  which  can  be  approximated  by  a  [Te  (r^f )].  If  in  eq.(47) 
the  pressure  p  is  now  replaced  according  to  the  above  relation  of  eq.(53),  and  T  =  Tg  by  eq.(26), 
one  obtains 


1  /f 

=  1-21  •  10-'®  •  r  • 


Herein  I  must  be  taken  in  [A],  the  mass  flow  rate  rh  in  [kg/s]  and  the  arc  constrictor  radius  Vc  in 
[m]  in  order  to  obtain  (we^e)rAf  dimensionless.  The  factor  /*  is  related  to  /  by 


where  for  n  ^  4  the  relation  must  be  considered  as  an  approximation. 

/•  now  is  a  slowly  varying  function  of  T*  and  p  or  of  ^  and  p,  or  of  ^  and  rh  according  to  eq.(53). 
In  Fig.20  the  factor  /*  is  plotted  for  argon  (/jj),  as  a  function  of  ^  with  the  pressure  as  a  parameter 
for  n=4.  One  can  see  that  within  the  ranges  of  interest  (  2.6  •  10*^  <  ^  <  4  •  10*^  ;  lO^Po  < 
p  <  10^ Pa  )  for  both  propellancs  /*  varies  only  very  little  and  that  within  these  ranges  /*  can  be 
approximated  by  a  constant  of  about  0.40  with  an  error  of  about  ±  0.06  for  argon.  Hence  it  is  a  fair 
approximation  to  replace  (/*)  within  eq.(54)  by  0.40.  The  product  of  average  cyclotron  frequency 
and  electron  collision  time  at  the  radius  rjvr  of  the  discharge  channel  where  the  self  magnetic  field 
becomes  a  maximum  follows  therefore  by 


1  li 

~  0.50  •  10-*®>/Aio  •  re‘  -A  (56) 

This  approximate  formula  shows  that  the  number  primarily  proportional  to  the  |  -th 

power  of  the  current  and  inversely  proportional  to  the  mass  flow  through  the  discharge  channel. 
With  increasing  chsmnel  or  constrictor  radius  r^,  {ytTe)rM  should  also  slowly  increase;  however, 
this  effect  is  very  weak  since  the  Tc  -dependency  is  only  proportional  to  the  5-th  power.  The 
atomic  weight  dependency  is  proportional  to  the  square  root  since  within  the  discharge  channel  the 
molecules  are  dissociated.  For  augon  as  propeUant  Mo  =  40. 
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Figure  20:  The  factor  /’  =  for  argon  as  a  function  of  -j-  with  p  as  a  parameter 


3.3  Experiment  and  Comparison  with  Calculations 

As  mentioned  above  the  magneto-plasma-dynamic  thruster  reaches  an  unstable  mode  at  high  power 
levels.  The  experiments  show  that  voltage  versus  current  curve  rises  steebly  above  a  certmn  current 
and  that  somewhat  above  this  threshold  voltage  fluctuations  occur  the  amplitudes  of  which  increase 
dramatically  until  the  thruster  becomes  irreversibly  damage.  The  ‘onset  current’  {Icrit)  is  now  taken 
as  that  current  at  which  the  voltage  noise  becomes  noticeable  but  at  which  the  thruster  is  still  not 
irreversibly  damaged.  The  measured  onset  conditions  so  far  are  plotted  in  Table  2.  Regarding  the 
DT-2  thruster  with  a  constrictor  diameter  of  24mm  the  onset  conditions  are  covered  between  a 
mass  flow  rate  of  0.3  g/s  and  1.5  g/s.  In  the  case  of  the  DT-5-thruster  with  a  constrictor  diameter 
of  30mm,  only  two  currents  at  a  mass  flow  rate  of  0.8  and  1.6  g/s  have  been  measured  with  argon 
as  propellant. 

For  the  DT-6  thruster  which  had  a  cylindrical  shape  of  constant  cross  sectional  area  for  the  arc 
chamber  and  the  throat  with  a  diameter  of  36mm,  the  onset  regime  is  couvered  between  a  mass 
flow  range  of  0.8  g/s  to  2.4  g/s  argon.  In  the  last  experiment  a  current  of  8300  A  has  been  reached 
at  only  a  faint  voltage  noise.  In  order  to  avmd  thrusterdamage  the  current  has  not  been  increased 
further. 

All  these  measurements  are  subtte  since  only  a  minor  current  increase  beyond  this  onset  can  lead 
to  a  faUore  and  total  damage  of  the  thruster.  E^pedally  at  higher  critical  currents  and  mass  flower 
rates  one  has  to  be  extremely  carefull  to  avoid  local  overheating  of  the  anode. 
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l^ble  2:  List  of  Available  Onset  Conditions  for  Different  Thruster  Types,  with  argon  as 

propellant. 


More  recent  measurement  with  a  CCD-camera  which  allows  exposure  times  of  100  ns  indicate  that 
near  onset  conditions  the  discharge  cannot  be  considered  being  rotational  symmetric  filling  out  the 
entire  discharge  chamber  or  nozzle  of  the  thruster,  but  rather  assumes  a  more  confined  channel 
configuration  which  moves  around  rapidly.  On  increasing  the  current  only  by  a  small  amount  above 
onset  the  anode  attatchment  tends  to  stick  more  and  more  to  a  local  point  thereby  overheating 
and  irreversibly  dama^ng  the  anode  at  a  confined  spot  (Fig.  21). 


Figure  21:  CCD-eamera  {dcture  at  onset  condition 
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For  the  sake  of  comparison  of  the  experiments  with  the  analysis  for  all  three  thruster  types  (DT- 
2,  DT-5  and  DT-6)  the  onset  conditions  -  Icrit  ^  a  function  of  rh  -  following  from  eq.(54)  were 
calculated  according  to  the  requirement 


7,65- 10 


-10 


A  /? 

rric 


ertt 


(57) 


with  {u}eTt)crit  'Values  of  0.5  and  0.75.  The  results  of  these  calculations  are  plotted  as  lines  for  the 
DT-2,  DT-5  and  DT-6  thrusters  in  Fig.  22,  Fig.  23  and  Fig.  24  respectively.  Herein  the  factor  /JJ 
was  determined  for  n=4  by  accounting  for  the  weak  dependency  of  m  and  One  can  see  that  the 
measured  points  lie  below  the  upper  curve  with  ~  ^  could  be  expected  and  that  the 

functional  dependency  of  the  onset  condition  for  argon,  ^ven  by 


2 

r|  =  1.95  •  10® 

me 


mii4s 


kg/s 


(58) 


are  within  the  accuracy  of  experiment  and  theory  (<  10%)  for  both  thruster  types.  According  to 
this  half  empirical  relation  (eq.(58))  should  be  set  equal  to  0.62  for  all  three  thrusters. 
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Figure  24:  Critical  current  as  a  function  of  the  mass  flow  rate  (DT-6) 
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3.4  Stability  Considerations 


The  previous  sections  have  shown  that  ‘onset’  can  be  quantitatively  explained  by  a  dramatic  decay 
of  the  radial  heat  conduction  loss  of  the  electrons  caused  by  the  azimuthal  self  magnetic  field.  When 
the  product  of  cyclotron  frequency  and  average  collision  time  of  the  electrons  exceeds  a  critical 
value  i{ueTe)rM  ^  0-®)  ^  noticeable  quenching  of  the  heat  conduction  occurs  leading  to  an  electron 
temperature  increase  within  the  center  of  the  discharge  channel  and  therefore  to  a  constriction  of  the 
electrically  conducting  channel.  This  run-away  heating  process  in  the  discharge  core  proceeds  until 
a  new,  energetically  balanced  discharge  channel  will  have  adjusted  itself.  This  new  adjustment 
is  characterized  by  the  fact  that  now  the  Joule  heat  is  balanced  by  the  heat  conduction  loss  of 
the  heavy  particles.  Since  the  conduction  coefficient  due  to  the  heavy  particles  is  smaller  by  a 
factor  of  y/me/mo  than  that  of  the  electrons  without  UeTe—  effect,  the  radial  temperature  slope 
should  steepen  and  consequently  lead  to  a  temperature  increase  in  the  discharge  channel  core.  This 
concentration  of  the  discharge  channel  now  causes  a  portion  of  the  gaseous  propellant  to  stream 
around  the  core  and  not  through  the  actual  current  carrying  channel.  The  mass  flow  rhc  flowing 
through  the  discharge  therefore  decreases  when  the  WeTe-  effect  sets  in  [{{LJeTe)rM  ^  0.6)]. 

Now  it  shall  be  shown  that  this  transition  from  a  high  to  a  low  heat  conduction  or  from  a  larger 
to  a  lower  mass  flow  rhc  and/or  from  a  chamber  filling  to  a  more  constricted  discharge  channel 
initiates  the  proper  instability.  So  far  only  a  straight  rotational  symmetric  discharge  configuration 
within  the  constrictor  has  been  considered.  However,  when  the  onset  occurs,  a  constricted  discharge 
channel  will  also  continue  dowr'  .  r  jn  into  the  diverging  nozzle  section  where  rotational  symmetry 
can  no  longer  be  expected.  Now  in  order  to  determine  the  new  ‘onset’  channel  configuration  from 
the  constrictor  downwards  to  the  anode  ring  a  previously  developed  theory  on  arc  stability  [29,  32] 
will  be  applied.  According  to  this  theory  the  configuration  of  a  current  carrying  plasma  channel 
is  only  a  stable  one  if  a)  the  sum  of  all  transverse  acting  forces  per  unit  channel  length  (xx) 
is  zero  along  the  entire  length  of  the  discharge  axis  and  b)  if  for  any  devitation  of  such  a 
balanced  configuration  a  resultant  force  Xx  /  0  is  initiated  pointing  in  the  opposite  direction  of 
the  deviation,  i.e.  tends  to  bring  the  channel  into  its  original  balanced  state. 

For  the  problem  considered  here,  this  force  accounting  for  all  body-  and  surface  forces  (e.g.  magnetic 
forces  caused  by  the  self  magnetic  field,  friction,  convection  and  pressure  effects  and  if  necessary 
Coriolis  and  centrifugal  forces)  follows  as  [29] 

Xx  =  [1  +  /  (n)]  -  pvldAc  -1-  pUzdAcj  ^  (59) 

Now,  R  in  eq.(59)  is  the  local  radius  of  curvature  (see  Fig.  25),  R  the  vector  pointing  to  the  center  of 
curvature  in  the  x-direction,  /  (n)  is  a  slowly  varying  function  of  the  grade  n  and  is  approximately 
equal  to  0.8  and  Vcy  is  the  transverse  velocity  of  the  channel  in  the  y-direction,  i.e.  normal  to  the 
axis  (z-direction)  and  normal  to  R  (see  Fig.  25). 

The  quantity  ^  is  the  so  called  “twist”  of  the  channel  and  follows  for  the  most  general  channel 
segment  configuration,  the  conical  helix  (see  Fig.  26)  by 

=  — sin/J  cosfiJl  —  tari^a  tan^0  (60) 

as  Th  ^ 
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Figure  26:  Right-tum  Conical  Helis  of  Constant  /?  and  a 

For  such  a  channel  segment  configuration  the  radius  of  currature,  R,  is  now  related  to  the  helix 
radius  and  the  angles  a  and  by 


(61) 
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Now,  the  axial  mass  transport  given  by 


/  pv^dAc  =  ±mc  (62) 

Jac 

is  positive  if  the  mass  flow  along  the  channel  axis  is  in  the  same  direction  as  the  electric  current 
and  negative  if  in  the  oposite  direction. 

According  to  geometrical  considerations  the  transverse  channel  velocity,  Vcy,  can  be  related  to  the 
rotational  frequency,  a;/,,  of  the  helix.  It  is 

Vcy  =  -rh.oJhsin^\J  1  +  tan^a  tan'^/3  (63) 

Herein  can  be  positive  or  negative  depending  on  wether  it  is  a  right  turn  revolution  with  respect 
to  the  main  helixal  axis,  z,  or  a  left  turn  revolution  with  respect  to  z  (seeFig.26).  By  means  of  the 
equations  (60),  (61),  (62)  and  (63)  the  formula  for  the  resultant  transverse  force,  eq.(59),  can  now 
be  written  in  the  form 


Xi  =  [1  +  /  (")]  ~  pvldAc  -  (±mc)  ThUhsin^  ^  (64) 

If  one  now  a.s.sumes  that  the  axial  speed  Vz  within  the  Joule  heated  channel  reaches  the  speed  of 
sound  (see  also  eq.(53)]  the  axial  impuls  transport  can  be  approximated  by 


(65) 


or  by  replacing  T  =  Tg  according  to  eq.(26)  one  obtains  after  some  rearangements 


Xi  = 


r,  .  g.  ^  1,95-10^°  hMl+a) 

l  +  /  n)  me  —  <  —ri  -  ,  ^  ,  W - rj - 

XtcJ  1  me  1  +  /  (n)  V  Mo 

r  rhUhsin/S  ianisj  1  +  tan^a  cos^l3 

1  +  f{n)  \  I  J  ^ 


Herein  I  must  be  taken  in  [A];  rc,R  and  r^^  in  [m],  Uh  in  [s”‘]  and  the  mass  flow  rate  through 
the  current  carrying  channel,  rhe  in  [feya”*]  in  order  to  obtain  the  resultant  force  per  unit  channel 
lenght  in  [iVm~*] .  Depending  if  the  mass  flow  and  the  obtain  current  have  the  same  or  the  opposite 
orientation  one  had  to  take  the  minus  of  the  positive  sign  respectivly.  If  one  applies  this  equation 
to  the  discharge  flow  of  our  MPD-thrusters,  the  electric  current  is  opposite  to  the  flow  (if  the  anode 
is  located  downstream);  therefore  the  plus  sign  in  eq.(66)  is  valid.  On  the  other  hand  0  must 
fall  in  the  range  -x  <  0  <  0{oTir  <  0  <  2ic,3k  <  0  <  47r,etc.)  such  that  sm/3  <  0.  Taking  a 
current-density  paraboloid  of  grade  n=4  and  argon  as  propellant  with  an  atomic  weight  Mo  =  40, 
the  transverse  force  per  unit  channel  lenght  becomes  proportional  to 
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Xx  ~  I  ~  7  •  lO^yJ'^A  (1  +  a)  -1,1-  lO^r^cJfc  (t)  tanfi\J\  +  tan^ot  cos'^^ 

(67) 

Herein  ianl3  is  larger  zero  <  /3  <  -^  ;  tt  <  /3  <  |7r ;  3ir  <  /3  <  ^■K,etc^  for  a  right  turn  helix 

and  negative  <  (3  <  0;  \-k  <  ^  <  2x;e<cj  for  a  left  turn  helix.  By  means  of  this  vector  ex¬ 
pression  the  plasma  behavior  within  the  given  MPD-thrusters  (see  Fig.(  ))  under  onset  conditions 
are  discussed. 


3.5  Plasmabehavior  under  Onset  Conditions 


As  long  as  the  expression  stays  well  below  1,95  •  10®  •>  the  UeTg-  effect  on  the 

conduction  heat  flux  of  the  electrongas  is  still  negligible  for  an  argon  MPD-flow.  In  this  case 
the  electrically  conducting  plasma  is  assumed  to  fill  up  the  volume  of  the  throat  (on  cylindrical 
part  of  the  fluxed  thruster)  and  the  diverging  part  of  the  thruster  nozzle.  The  discharge  axis 
should  therefore  cioncide  with  the  center  line  of  the  thruster  and  any  deviation  of  the  straight 
configuration  should  be  meaningless  since  the  discharge  channel  is  “anchored”  at  the  cathode  tip 

and  fixed  by  the  chamber  walls.  If,  however,  comes  close  or  even  exceeds  the  critical 

[I  1 1 

run  away  heating  and  a  concentration  of  the  discharge  channel  sets  in. 
Simultaneously  with  the  concentration  one  must  expect  an  increase  in  the  characteristics  quantity 
^  since  now  a  larger  part  of  the  gasflow  streams  outside  of  the  hot  more  confined  discharge 

channel  yet;  this  part  cannot  be  accounted  any  more  to  rhc.  As  a  consequence  the  quantity  ^ 

r  9  i  1 

may  increase  even  beyond  a  value  of  1 , 7- (1  +  «)  >  2, 3-10®  what  caused  a  dramatic 

change  in  the  direction  of  the  resultant  transverse  force,  xx  which  acts  on  the  discharge  channel  per 

unit  lenght  [see  eq.(67)].  Based  on  this  proportionality  formula,  eq.(67),  for  xx  the  following  three 

(a  i\  r  a  1' 

1  ^  1,95-10®  'J-  . 


m 


<  i,7-io®vWiT^ 


ib)  =  1,7 -10® 

C*")  (^r|)<  1,7- 10®  v^MIT^O 


In  case  (a)  the  force  xx  points  into  the  negative  R—  direction.  If  an  anchorde  discharge  channel 
experiences  any  small  deviation  (J  ^  see  Fig.(27),  the  force  xx  tends  to  straighten  the  column. 
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In  this  case,  therefore,  a  straight  discharge  channel  is  stable.  Within  the  diverging  part  of  the  nozzle 
the  overheated  constricted  channel  must  be  curved,  however,  in  order  to  reach  the  anode  ring.  The 
raxlius  of  curvature  becomes  finite,  and  in  order  to  fulfill  the  requirement,  xj.  =  the  brace 
expression  of  eq.(67)  must  be  zero  or  the  channel  must  be  twisted  i.e.  assumes  within  the  diverging 
part  of  the  nozzle  (see  Fig.28). 

disturbed 


Cathode  straight  channel 


Figure  27:  Dlustration  of  a  straight  and  disturbed  anchored  discharge  channel 

In  thi'  case  for  a  right  turn  helix  {tana  >  0)  must  be  negative  or  for  a  left  turn  helix  {tana  <0)uk> 
0  in  order  to  achive  that  the  transverse  force  xx  becomes  zero. 

If,  however,  case  (b)  is  reached,  no  twisted  channel  is  necessary  in  order  to  obtain  Xx  =  0  and 
uiii  can  be  zero.  This  means  that  the  anode  attachement  sticks  at  a  point  and  damages  the  anode 
locally  by  overheating  which  can  indeed  be  observed  when  the  preonset  conditions  are  exceeded. 

In  case  (c)  the  transverse  force  Xx  per  unit  channel  lenght  (see  eq.(67))  is  zero  if  the  channel  is 
straight  {R  — >  oo).  Any  small  deviation  of  the  anchored  arc  channel  (see  Fig.27),  however,  now 
leads  to  an  unstable  channel  behaviour  or  eventually  to  a  curved  and  twisted  channel.  According 
to  eq.(67)  for  a  twist-free  but  disturbed  (R  finite)  channel  the  brace  expression  is  greater  than  zero, 
hence  Xx  acts  in  the  direction  of  R  ^  and  tends  to  disturb  the  channel  even  further.  Again  a 
belanced  configuration  with  Xx  =  0  can  only  be  archived  if  the  brace  expression  becomes  zero,  i.e. 
if  the  third  term  within  the  brace  of  eq.(67)  becomes  negative.  In  that  case  the  discharge  channel 
assumes  a  rotating  helical  shape  whereby  Uh  must  be  larger  than  zero  for  a  right  turn  helix,  and 

9  ^ 

than  zero  for  a  left  turn  helix.  It  is  interesting  to  note  here  that  by  increasing  the  quantity  the 
helix  should  again  straighten  itself  since  now  the  brace  expression  becomes  positive  and  therefore 
Xx  acts  in  the  direction  of  R  and/or  the  brace  expression  becomes  zero  which  means  tan^  increases 
and  therefore  brings  the  discharge  channel  to  a  straighter  configuration.  Within  the  diverging  part 
of  the  nozzle  the  helix  continues  in  a  conical  fashion  in  order  to  end  up  on  the  anode  ring.  Since 
must  be  unequal  zero,  the  helix  and  the  anode  attachment  rotates  again,  now  however,  in  the 
opposite  direction  than  in  case  (a). 
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(a)  “side— on"  (b)  “end-on“ 

Figure  28:  Helical  discharge  channel  within  nozzle 


Summarizing  it  can  be  said  that  at  ‘onset 


,  (M 

\  me 


10®  f 


L 


Aims 

kg/a 


the  discharge  concentrates 
This  conical  helix  rotates 


and  assumes  a  helical  shape  within  the  diverging  part  of  the  nozzle 
as  long  as  the  critical  condition  is  not  exceeded.  If  now  the  current  is  increased  and/or  the  mass 
flow  decreased  beyond  onset  conditions,  the  rotation  ceases  and  the  anode  attachment  begins  to 
stick,  thereby  damaging  the  anode  irreversibly  by  local  overheating.  By  further  increasing  the 

9  ^ 

characteristic  quantity  ,  the  helix  starts  to  turn  again  (now  in  the  opposite  direction  than 
before).  Since  under  these  extreme  conditions  a  straight  discharge  channel  is  unstable,  plasma 
fluctuations  near  the  cathode  and  within  the  constrictor  should  arise  which  may  now  damage  the 
inner  thruster  walls. 
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